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1.  INTRODUCTION 


A  steady-state  solution  q  of  the  two-dimensional  incompressible  continuity  and  Navier- 
Stokes  equations  which  describe  flow  in  a  prescribed  two-dimensional  domain  fl  bounded 
by  dCt  is  seeked  numerically.  A  plethora  of  numerical  approaches  for  the  accurate  and 
efficient  integration  of  either  the  steady  or  the  unsteady  equations  of  motion  exists  (e.g.  [9, 
26, 30])  so  that  this  problem  may  be  considered  solved  in  principle.  However,  in  performing 
a  time-accurate  integration  of  the  equations  of  motion  one  observes  that,  depending  on  the 
values  of  parameters  such  as  the  flow  Reynolds  number,  in  the  limit  of  large  time  either 
a  steady-state  solution  is  obtained  (e.g.  [1])  or  unsteady,  sometimes  periodic,  motion  sets 
in  (e.g.  [13,  11]).  The  first  question  arising  is  what  type  of  physical  information  is  not 
considered  by  solving  the  steady  as  opposed  to  the  unsteady  equations  of  motion  and  what  is 
the  physical  interpretation  of  the  critical  conditions  beyond  which  the  steady  and  unsteady 
formulations  deliver  different  results.  Both  physical  and  numerical  experience  suggest  that 
at  low  Reynolds  numbers  the  two  formulations  may  be  used  interchangeably.  For  example, 
essentially  identical  results  with  those  of  Ghia  et  al  [27]  and  Schreiber  and  Keller  [19] 
have  been  obtained  by  a  multitude  of  subsequent  investigators  who  used  the  time-dependent 
equations  of  motion  to  describe  flow  in  the  square  lid-driven  cavity  at  Reynolds  numbers 
up  to  Re  =  104.  On  the  other  hand,  the  question  of  existence  of  a  steady-state  solution 
delivered  by  the  unsteady  version  of  the  equations  of  motion  at  Re  —  104  has  been  recently 
re-opened  [31],  while  it  is  well  known  that  Hopf  bifurcations  exist  in  both  the  aspect-ratio 
two  singular  lid-driven  cavity  at  Re  <  5000  [11]  and  its  regularised  square  counterpart  at 
Re  «  104  [10].  Concensus  exists  that  at  high  Reynolds  numbers  the  unsteady  formulation 
is  capable  of  delivering  physics  inaccessible  to  the  steady  version  of  the  equations  of 
motion;  however,  the  origin  of  the  differences  between  the  results  of  the  two  formulations 
is  presently  not  understood  in  a  satisfactory  manner.  This  is  an  alternative  way  of  posing 


the  first  question,  namely  what  are  the  unsteady  effects  that  manifest  themselves  at  high 
Reynolds  numbers. 

The  next  question  arises  from  the  very  concept  of  two-dimensionality.  The  results  of 
numerical  solutions  of  the  three-dimensional  analoga  of  the  incompressible  continuity  and 
Navier-Stokes  equations  are  in  most  cases  in  substantial  qualitative  and  quantitative  dis¬ 
agreement  with  their  two-dimensional  counterparts  (e.g.  [2, 7]),  relegating  two-dimensional 
DNS  to  the  realm  of  academic  interest.  Within  the  scope  of  two-dimensional  solutions  be¬ 
ing  of  interest,  three-dimensionality  of  physical  space  could  be  addressed  by  considering 
the  flow  to  be  independent  of  the  third  spatial  direction.  Homogeneity  in  this  third  direction 
could,  in  turn,  be  discussed  in  the  context  of  a  three-dimensional  simulation,  nonperiodic 
in  the  same  two  spatial  directions  as  the  two-dimensional  one  and  periodic  in  the  third. 
Advances  in  both  algorithms  and  hardware  and,  not  least,  a  considerable  amount  of  knowl¬ 
edge  on  the  differences  between  two-  and  three-dimensional  numerical  simulation  results 
lead  one  to  employ  a  DNS  algorithm  for  flow  with  two  nonperiodic  and  one  periodic  spatial 
direction  (e.g.  [22,  15,  28])  in  the  founded  expectation  that  a  three-dimensional  so-called 
spatial  DNS  is  the  only  means  capable  of  capturing  all  physical  phenomena  at  a  certain 
Reynolds  number.  The  second  question  which  may  be  posed  at  this  point  regards  the  origin 
of  the  differences  between  the  results  of  such  two-  and  three-dimensional  direct  numerical 
simulations.  Associated,  one  may  ask  whether  there  exists  an  alternative  means  to  spatial 
DNS  for  the  description  of  the  origins  of  the  three-dimensional  phenomena  encountered. 

The  objective  of  the  present  paper  is  to  put  both  questions  within  the  unified  framework 
of  nonparallel  linear  instability  of  the  steady  state  q.  With  the  aid  of  a  well-studied  flow 
example  we  demonstrate  the  intimate  link  between  numerical  residuals  in  steady-state  fluid 
flow  calculations  and  linear  two-dimensional  eigenmodes  of  the  converged  steady  state 
q.  In  §2  we  present  theoretical  arguments,  first  analysing  the  behaviour  of  numerical 


residuals  near  convergence  towards  the  steady-state  solution  from  a  numerical  point  of 
view.  Subsequently  we  discuss  solutions  of  the  partial  derivative  eigenvalue  problem 
governing  linear  instability  of  nonparallel  two-dimensional  steady-state  flows,  which  shed 
light  on  residuals  from  a  physical  viewpoint.  With  the  origin  of  residuals  established  from 
a  physical  point  of  view  we  construct  and  present  an  algorithm  which  permits  recovery  of 
the  converged  steady-state  solution  from  transient  results  of  the  time-marching  procedure, 
the  latter  taken  well  before  convergence.  In  §3  we  present  results  obtained  for  the  classic 
square  lid-driven  cavity  flow  as  a  demonstrator  of  the  ideas  discussed  herein.  The  link 
between  the  results  of  nonparallel  linear  instability  theory  and  different  types  of  behaviour 
of  numerical  residuals  in  the  DNS  is  demonstrated  in  this  section  and  the  aforementioned 
questions  are  answered.  Examples  of  recovery  of  the  converged  steady-state  from  transient 
data  and  assessment  of  the  substantial  savings  in  the  computing  effort  materialised  by  use 
of  the  proposed  algorithm  are  presented  in  this  section.  Closing  remarks  on  the  far-reaching 
implications  of  the  present  findings  are  made  in  §4  and  suggestions  for  the  extension  of  the 
present  analysis  to  compressible  flow  and  flow  with  three  nonperiodic  spatial  directions  are 
made  in  §5. 


2.  THEORY 

2*1’  On  residuals  and  the  phenomenology  of  their  behaviour 

While  in  an  computation  based  on  the  steady  system  of  equations  governing  fluid  flow 
motion  residuals  are  viewed  as  departure  from  the  steady  state  which  have  to  be  eliminated 
in  an  efficient  manner  by  a  specific  solution  algorithm  (e.g.  multigrid),  in  a  time-accurate 
integration  one  may  view  transients  as  solutions  of  the  equations  of  motion  and  attempt 
to  attach  physical  significance  to  characteristic  patterns  of  their  behaviour.  Here  we  con¬ 
centrate  on  a  time-accurate  integration  of  the  unsteady  equations  of  motion  and  monitor 
the  behaviour  of  residuals,  defined  as  the  difference  between  the  transient  solution  and  the 


converged  steady  state,  in  flow  regimes  where  the  latter  exists.  Physical  space  is  three- 
dimensional;  without  loss  of  generality  we  may  take  the  Cartesian  coordinates  x  and  y  to 
be  defined  on  H  while  z  denotes  the  third  spatial  coordinate  in  the  direction  of  ft.  Along 
the  first  two  coordinates  the  velocity  vector  has  components  u  and  v,  while  pressure  is 
denoted  by  p.  The  equations  of  motion  are  marched  in  time  t  until  q  =  (u,v,p)T,  the 
transient  solution,  converges  to  q.  Assuming  that  the  latter  exists  and  keeping  the  domain 
ft  unchanged,  the  following  qualitative  observations  are  made. 

First,  at  any  Reynolds  number  Re  at  which  q  exists,  close  to  convergence  the  residuals 
decay  exponentially  in  amplitude.  Second,  refinement  of  the  discretisation  of  the  domain 
ft  at  constant  Re  results  in  convergence  of  the  rate  at  which  the  residuals  decay.  Third,  the 
(converged)  rate  of  decay  of  residuals  is  a  function  of  the  flow  Reynolds  number;  as  Re 
increases  residuals  decay  slower  and  the  associated  time  of  integration  of  the  equations  of 
motion  until  convergence  increases.  Fourth,  on  occasion,  the  residuals  decay  at  a  specific 
constant  rate  for  a  number  of  decades  before  this  rate  of  decay  changes  to  a  different 
constant  value  at  which  residuals  further  decay  until  convergence.  Fifth,  systematically 
increasing  Re ,  instead  of  monotonic  convergence  of  residuals  an  oscillatory  behaviour  of  q 
in  the  neighbourhood  of  q  is  observed.  Ultimately,  a  value  of  Reynolds  number  is  reached 
past  which  no  q  exists.  At  first  sight  the  existence  of  a  physical  mechanism  which  unifies 
such  diverse  patterns  of  behaviour  of  the  numerical  solution  seems  unlikely. 

2.2.  A  numerical  point  of  view  on  the  behaviour  of  residuals  near  convergence 

However,  it  is  straightforward  to  provide  an  explanation  of  the  first  observation  on  the 
behaviour  of  residuals,  which  also  provides  a  handle  to  the  link  between  numerical  residuals 
and  physical  flow  instabilities.  We  assume  that  the  solution  q  is  close  to  converging  to 
the  seeked  two-dimensional  field  q  =  (u,fi,p)T  such  that  it  may  be  decomposed  into  the 
latter  and  small  two-dimensional  residuals  q2D  =  (^2D,  V2d,P2d)t  superimposed  upon  it, 


according  to 


q(x,y,t)  =  q(s,j/) +eq2D(*,y,*).  (!) 

with  £«1,  We  next  substitute  the  decomposition  (1)  into  the  continuity  and  Navier- 
Stokes  equations  and  assume  that  the  steady-state  solution  satisfies  the  equations  of  motion 
at  0(1) ,  such  that  it  may  be  subtracted  out  of  the  resulting  system  at  this  order.  Subsequently, 
based  on  the  smallness  of  the  amplitude  of  the  residuals,  we  linearise  about  q  and  rearrange 
the  system  at  0(e)  such  that  the  vector  of  residuals  represents  the  unknowns;  terms  of 
0(e2)  are  neglected.  Since  the  coefficients  of  the  resulting  linear  system  of  equations  for 
the  determination  of  q2D  at  0(e)  are  independent  of  time  t  we  may  introduce  an  eigenmode 
decomposition  in  this  coordinate,  according  to 

q2D(*,y,<)  =  q2D(z,y)  e0-4  (2) 

with  q2D  =  («2Dj #2D,P2d)t-  The  physical  significance  of  the  parameter  a  will  be 
discussed  shortly;  from  a  numerical  point  of  view  it  represents  the  rate  at  which  the 
residuals  q2D  decay  in  the  neighbourhood  of  q.  For  simplicity  we  present  only  the  real  part 
of  the  admissible  solutions  of  (2)  although  it  is  clear  that  both  q2p  and  a  may,  in  general, 
be  complex  while  q2D  is  always  real.  Convergence  of  the  solution  q  towards  q  may  be 
monitored  by  reference  to  either  the  local  behaviour  of  the  solution  q  at  a  position  (z0 ,  yo) 
on  ft  or  by  monitoring  a  suitably  defined  global  criterion  such  as  the  energy  contained 
in  the  residuals  q2n;  alternatives  have  been  discussed  in  [24].  Here  we  follow  the  first 
approach  and  recover  the  parameter  a  by  monitoring  the  solution  at  two  time-levels,  t  —  At 
and  t ,  where  At  may  but  need  not  be  the  time-step  in  the  numerical  solution  algorithm. 
Combining  (1)  and  (2)  it  follows  that  the  time-behaviour  of  the  solution  may  be  monitored 
by 


cr  =  ln[q</q<  At)/At  «  dlnfq^/dt, 


(3) 


where 


q‘ =  |q(x0lyo,<)  -  q(xo,yo)|-  (4) 

The  approximation  in  (3)  holds  as  equality  in  the  case  of  linear  dependence  of  lnfq*]  on 
time  t.  Decay  of  residuals  is  indicated  by  a  <  0.  The  first  statement  of  the  present  paper 
is  thus  in  place  without  reference  to  a  particular  flow,  through  the  analytical  result  that  an 
exponential  decay  of  residuals  near  convergence  should  be  observed  as  a  consequence  of 
the  separability  of  the  linearised  system  of  equations  for  the  determination  of  residuals  in 
time. 

2.3.  A  physical  point  of  view  based  on  nonparallel  linear  instability  theory 

Explanation  of  the  further  observations  made  in  §  2.1  requires  calling  upon  an  extension 
of  the  classic  linear  instability  theory  proposed  by  Tollmien  [25],  which  describes  the 
behaviour  of  small-amplitude  disturbances  superimposed  upon  an  one-dimensional  steady- 
state  basic  profile,  into  a  new  theory  which  is  concerned  with  small-amplitude  perturbations 
superimposed  upon  a  steady  two-dimensional  field.  In  so  doing,  the  many  and  often 
questionable  assumptions  related  with  the  so-called  parallel-flow  approximation  are  relaxed 
and  the  linear  instability  of  nonparallel  basic  states  may  be  analysed.  The  penalty  to  be 
paid  in  resolving  two  spatial  dimensions  is  the  need  for  numerical  solution  of  a  partial- 
derivative-based  eigenvalue  problem  instead  of  the  straightforward  ordinary-differential- 
equation-based  system  of  the  Orr-Sommerfeld  and  Squire  equations  [6].  One  of  the  early 
successes  of  the  nonparallel  two-dimensional  linear  instability  analysis  was  the  discovery 
of  inviscid  short-wave  instability  of  two-dimensional  eddies  by  Pierrehumbert  [18]  while 
the  first  viscous  linear  analysis  in  two  non-periodic  spatial  dimensions  known  to  us  is  the 
work  of  Lee  et  al.  [16]  on  the  instability  of  flow  in  a  rectangular  enclosure  under  the 
influence  of  gravity  and  temperature  gradient.  More  recent  viscous  analyses,  in  step  with 


modem  developments  in  algorithms  and  hardware,  have  been  presented  in  [23]  where  we 
discuss  the  linear  instability  of  two  nonparallel  flows,  that  in  a  rectangular  duct  and  that 
in  the  infinite-swept  attachment-line  boundary  layer.  The  nonparallel  linear  instability  of 
laminar  flows  encompassing  a  recirculation  bubble  has  been  discussed  by  Barkley  et  al 
[4]  for  flow  behind  a  backward-facing  step  and  by  Theofilis  et  al  [29]  for  two  separated 
boundary  layer  flows  developing  under  the  influence  of  adverse  pressure  gradients. 

We  re-interpret  the  transient  solution  q  in  three-dimensional  physical  space  as  one  com¬ 
posed  of  small-amplitude  three-dimensional  perturbations  q  =  (u,v,w)p)T  superimposed 
upon  q  =  (u,  v ,  tu,p)T,  the  latter  again  taken  to  be  two-dimensional.  Linearisation  about 
q  is  permissible  on  account  of  the  smallness  of  perturbations  compared  with  the  steady- 
state  q  and  the  resulting  system  for  the  determination  of  q  is  separable  in  both  t  and  z  on 
account  of  the  steadiness  and  the  two-dimensionality  of  the  basic  flow  q.  Eigenmodes  are 
introduced  in  these  directions  such  that 

q(z,  y)  z,  t)  =  q(x,  y)  e1  +  c.c.  (5) 

with  q  =  (u,  v,  w,  p)T  and  w  being  the  disturbance  velocity  component  in  the  z-direction. 
Complex  conjugation  is  introduced  in  (5)  since  q  is  real  while  all  three  of  q,/3  and  uj  may 
be  complex.  In  the  framework  of  a  temporal  linear  nonparallel  instability  analysis  used 
presently  we  write  the  linearised  system  in  the  form  of  an  eigenvalue  problem  for  the 
complex  quantity  u,  while  0  is  taken  to  be  a  real  wavenumber  parameter  describing  an 
eigenmode  in  the  z-direction.  The  real  part  of  u>  is  related  with  the  frequency  of  the 
instability  mode  while  its  imaginary  part  is  the  growth/damping  rate;  a  positive  value  of 
Wi  =  $s{u;}  indicates  exponential  growth  of  the  instability  mode  q  in  time  t  while  <  0 
denotes  decay  of  q  in  time.  In  the  present  framework  the  three-dimensional  space  comprises 
0  extended  periodically  in  z  and  characterised  by  a  wavelength  Lz  in  this  direction  which 
is  associated  with  the  wavenumber  of  each  eigenmode,  0,  through  Lz  —2^/0. 


The  system  for  the  determination  of  lj  and  q  takes  the  form  of  a  complex  nonsymmetric 
generalised  eigenvalue  problem 


[C  -  (Vxu)\  it  -  ( T>yu)v  -  Vxp 

~(DXV)U  +  [C  -  {VyV)\  V  -  VyP 
~(Vxw)u  —  (T>yw)v  +  Cw  —  i  f3p 
Vxu  +  VyV  +  i  f3w 


—  -i  wu,  (6a) 

=  -i  u;  v,  (6b) 

—  -i  u  u),  (6c) 

=  0,  (6d) 


subject  to  appropriate  boundary  conditions  on  50.  The  linear  operator  C  =  (1  /Re) 
(V2X  +T>1~  (32)  -  uVx  -vVy-i  pw  with  Vx  =  e/dx,  Vl  =  d2 /dx2,  Vy  =  d/dy 
and  =  d2/ dy2 .  Some  details  of  an  efficient  numerical  algorithm  for  the  solution  of  the 
matrix  eigenvalue  problem  resulting  from  numerical  discretisation  of  the  spatial  directions 
x  and  y  have  been  presented  in  [23]. 

Comparison  of  (1-2)  and  (5)  reveals  that  the  two  formalisms  are  related  in  the  limit 
0  0.  However,  w  is  not  taken  a  priori  to  vanish  within  the  framework  of  nonparallel 

linear  instability;  three-dimensionality  of  physical  space  is  preserved  and  the  existence  of  a 
two-dimensional  steady-state  solution  q  is  the  result  of  q  ■  ->  0  as  t  — >  oo.  The  comparison 
of  (1-2)  and  (5)  highlights  two  further  key  ideas  of  the  present  paper.  On  the  one  hand, 
the  residuals  discussed  earlier  acquire  the  physical  interpretation  of  one  of  the  linear 
eigenmodes  which  pertain  to  the  steady-state  q  and  have  0  =  0;  on  the  other  hand,  the  rate 
of  decay  of  the  residuals  o  is  nothing  but  the  damping  rate  u-,  of  this  linear  perturbation, 
as  delivered  by  numerical  solution  of  the  partial-derivative  eigenvalue  problem  (6a-6d). 
Another  question  naturally  arising  concerns  the  physical  behaviour  of  the  system  when  the 
least  stable  member  of  the  linear  eigenspectrum  which  pertains  to  q  and  has  0  =  0  becomes 
unstable.  The  answer  is  clearly  that  the  existence  of  an  unstable  (0  =  0) -eigenmode  is 


mutually  exclusive  with  the  ability  to  obtain  a  converged  q.  From  the  point  of  view  of  the 
global  linear  instability  theory  based  on  the  partial  derivative  eigenvalue  problem  (6a-6d) 
the  unsteady  behaviour  of  two-dimensional  flow  may  be  related  to  (0  =  0)— eigenmodes 
approaching  conditions  of  neutral  stability  and  interacting  nonlinearly. 

The  answer  to  the  second  question  posed  in  the  Introduction  may  now  also  be  obtained 
without  reference  to  a  specific  flow  example.  The  existence  of  a  steady-state  q  in  a  2D 
numerical  simulation  is  synonymous  with  the  fact  that  all  (0  =  0) -eigenmodes  of  the  flow 
have  <  0.  Modes  having  0  /  0,  on  the  other  hand,  may  be  either  growing  or  decaying 
linearly.  In  case  <  0  V  0,  a  three-dimensional  numerical  simulation  performed  at  some 
parameters  in  a  three-dimensional  domain  defined  by  H  and  an  arbitrary  periodic  extent 
Lz  in  the  z-direction  will  deliver  identical  results  for  a  converged  q  compared  with  that 
of  a  two-dimensional  simulation  performed  at  the  same  parameters  in  the  domain  SI  The 
situation  changes  in  case  a  bracket  of  wavenumbers  0  £  [0i,02\  exists  which  corresponds 
to  unstable  modes.  The  largest  wavenumber  02  defines  a  length  Lz2  =  2tv/02\  if  the 
three-dimensional  simulation  is  performed  with  Lz  <  Lz2  again  no  difference  is  to  be 
expected  between  its  result  for  q  and  that  of  a  two-dimensional  simulation.  Both  will 
converge  to  the  same  steady-state  solution  q  since  all  wavenumbers  of  modes  defined  by 
an  Lz  constrained  as  above  correspond  to  u)x  <  0.  However,  if  Lz  >  Lz2  at  least  one 
mode  in  the  three-dimensional  simulation  will  be  unstable,  which  will  result  in  the  two- 
and  three-dimensional  simulations  producing  different  solutions. 

We  return  to  the  observation  of  oscillatory  behaviour  of  the  residuals  near  convergence 
and  differentiate  between  exponentially  decaying  residuals  of  either  sinusoidal  or  apparently 
nonlinear  nature.  A  linear  decay  of  In [q*]  is  a  consequence  of  (0  =  0) -linear  eigenmodes 
being  stationary,  i.e.  having  ujt  =  S£{u>}  =  0.  However,  other  stable  two-dimensional 
member  of  the  eigenspectrum  of  q  need  not  correspond  to  stationary  modes;  damped 


travelling  modes  having  vT  ^  0  will  manifest  themselves  in  the  time-accurate  simulation 
as  residuals  of  sinusoidal  character  the  magnitude  of  which  decays  exponentially.  On  the 
other  hand,  the  unambiguously  linear  dependence  of  lnfq*]  on  t  in  the  neighbourhood 
of  q  is  the  consequence  of  the  existence  of  a  spectrum  comprising  modes  which  are 
clearly  separated  in  parameter  space  from  one  another.  The  co-existence  of  several  two- 
dimensional  eigenmodes  of  approximately  the  same  damping  rate  can  lead  to  their  nonlinear 
interaction  and  difficulty  to  observe  a  behaviour  governed  by  nonparallel  linear  instability 
theory.  Comparison  of  power  spectral  analysis  of  the  time-dependent  DNS  signal  and 
the  results  of  the  partial-derivative  eigenvalue  problem  (6a-6d)  may  shed  light  upon  the 
two-dimensional  eigenmodes  involved  in  such  a  nonlinear  interaction. 

2.4.  On  the  time  of  integration  until  convergence 

Straightforward  rearrangement  of  (1-2)  delivers  an  estimate  of  the  time  necessary  (under 
linear  conditions)  for  the  least  stable  global  mode  present  in  the  numerical  solution  to  be 
reduced  from  an  amplitude  Ai  to  A2 ,  which  may  be  calculated  from 


TAi/A3=HAi/A0)/( —a*),  (7) 

where  is  the  damping  rate  of  the  mode  in  question.  The  worst  case  scenario  in  a  time- 
accurate  integration  is  that  the  solution  will  be  attracted  by  the  least-stable  global  eigenmode 
developing  upon  q  and  having  0  =  0  throughout  the  course  of  the  simulation.  An  upper 
bound  for  the  time  necessary  for  the  steady-state  to  be  obtained  may  then  be  offered  by 
(7)  in  which  ui  is  the  damping  rate  of  this  mode.  Defining,  for  example,  convergence 
as  the  reduction  of  an  0(1)  residual  by  10  orders  of  magnitude  results  in  an  integration 
time  of  T10-io  «  23/|u>i|.  This  is  a  conservative  estimate  since  it  is  occasionally  observed 
that  other  stronger  damped  eigenmodes  will  come  into  play  early  in  the  simulation  and  the 
least-damped  eigenmode  will  only  determine  the  late  stages  of  the  convergence  process. 


An  associated  point  concerns  the  misconception  which  often  exists  that  initialising  the 
numerical  solution  for  q  at  some  Reynolds  number  from  a  state  which  is  ’close’  to  the 
one  desired,  for  instance  using  the  converged  solution  at  a  somewhat  different  Reynolds 
number,  may  reduce  the  integration  time.  In  the  context  of  the  present  analysis  this  is  shown 
to  be  a  misplaced  expectation.  If  there  exists  an  0(1)  deviation  between  the  target  solution 
and  its  initial  estimate,  the  deviation  has  to  be  reduced  in  magnitude  during  an  integration 
of  the  equations  of  motion  for  the  length  of  time  determined  by  the  least  damped  two- 
dimensional  (p  =  0) -eigenmode  at  the  specific  Reynolds  number.  It  is  this  eigenmode  of 
the  flow  and  not  the  initial  state  which  determines  the  length  of  the  integration  time  for  q. 
The  ideas  discussed  in  §2.3,  on  the  other  hand,  lead  to  an  algorithm  application  of  which 
may  save  substantial  amounts  of  the  integration  time  necessary  for  reduction  of  residuals 
to  machine-roundoff  level. 

2.5.  Recovery  of  the  converged  solution  q  from  transient  data 

Having  identified  small-amplitude  residuals  in  the  calculation  as  the  least  damped  global 
two-dimensional  eigenmodes  of  the  flow,  it  is  now  possible  to  utilise  this  information  in 
order  to  recover  the  converged  steady-state  solution  from  transient  data,  without  having  to 
pursue  the  time  integration  of  the  equations  of  motion  until  convergence  in  time  is  obtained. 
Combining  (1),  (2)  and  (5)  one  obtains 

q(x,y,i)  =  q(z,y)  +  e  [q*  cos wTt  -  %  sin wTt j  efft,  (8) 

where  q,  =  &{q},  Qi  =  S{q}  and  q  is  one  of  the  (P  =  0) -linear  eigenmodes  in  (5). 
It  should  be  stressed  here  that  the  following  discussion  is  applicable  to  transient  data  for 
which  (8)  holds,  namely,  solutions  for  which  the  entire  time-dependence  of  the  solution 
is  exhibited  in  the  residuals;  in  other  words,  the  present  analysis  is  based  on  the  self- 
consistent  premises  that  3q jdt  =  0.  Further,  it  is  noted  that  q  may  but  need  not  be  the 


least-damped  member  of  the  eigenspectrum  of  q;  the  only  prerequisite  for  the  validity  of 
the  following  discussion  is  that  the  transient  solution  has  reached  a  regime  of  exponential 
decay  of  residuals.  A  final  point  is  that  the  signal  near  convergence  need  not  be  composed 
of  a  single  damped  eigenmode  as  (8)  implies.  However,  the  elements  of  the  theory  for  the 
recovery  of  q  from  a  signal  being  composed  of  several  stationary  (a;r  =  0)  and  travelling 
(a )T  ^  0)  linearly  damped  eigenmodes  may  be  exposed  by  reference  to  (8)  on  which  we 
focus  our  attention. 

The  calculation  of  q  from  transient  data  for  q  follows  in  two  stages.  First,  elementary 
signal  analysis  techniques  deliver  the  results  for  wT  and  a.  Second,  once  wr  and  a  have 
converged  in  time  (8)  may  be  used  to  calculate  q.  The  circular  frequency  wT  is  calculated 
from  the  the  period  of  oscillations  in  the  time-signal  of  q  which,  in  turn,  is  identified  by 
the  maxima  in  the  signal.  Independently,  in  order  to  calculate  a  we  re-write  (8)  as 


03q 
dt 3 


dt 


+  (a2 +^-2^=0. 


_d\ 

dt2 


(9) 


This  expression  may  be  evaluated  at  those  times  that  dq/ dt  =  0  in  the  course  of  the 
time-integration,  i.e.  at  the  same  times  that  a;r  is  calculated.  At  these  times  the  magnitude 
of  cr  is  given  by 


a 


1  (33q/S<3) 


2  (32q /dt2) 


(dq/dt)= 0 


(10) 


In  case  ujr  =  0,  a  monotonic  dependence  of  dq/dt  on  t  is  usually  observed  from  the 
beginning  of  the  calculation  until  convergence,  with  dq/dt  =  0  only  at  convergence.  In 
this  case,  the  magnitude  of  a  may  be  calculated  using 


(gya*2) 

(gq /dt) 


(11) 


With  cr  and  wr  converged  in  time  (8)  may  be  written  as  a  linear  system  of  three  equations 
at  three  times  h  =  t,  t2  =  t  +  At  and  t3  =  t  +  2A t  for  three  unknowns,  q,  qr  and  qi  with 


the  transient  solution  qn  =  q(x,  y,  tn)  known  at  these  times.  Simple  algebra  delivers  the 
desired  converged  steady-state  solution  q  as 

-  _  qi e2<rAt  z 2  q2  e<rAt  c0S(i;rAt  +  q3  ,l2) 

q  e2<rAt  _  2  COSU>rA£  +  1 

As  an  aside,  the  spatial  structure  (4r ,  4)  of  the  linear  eigenmode  q  may  also  be  recovered 
to  within  an  arbitrary  constant  from  the  same  linear  system.  Equivalently,  if  only  the 
converged  steady-state  solution  is  of  interest,  the  expression 


1  2,  «  Sq  S2q\ 


(13) 


may  be  used  for  the  recovery  of  q  from  transient  data  for  q  and  its  first  two  time- 
derivatives.  Either  of  (12)  or  (13)  may  be  used  for  the  cases  of  residuals  corresponding  to 
stationary  (o;r  =  0)  or  travelling  (u;r  ^  0)  single  linear  eigenmodes. 

This  idea  may  be  extended  to  extract  q  from  a  DNS  signal  comprising  several  linearly 
decaying  eigenmodes  superimposed  upon  the  steady-state  solution, 


q  =  q  +  ^£n(gn,r  cos c -  qn,i  sin  vrlyrt)e<Tnt .  (14) 

n 

As  an  example,  in  the  case  of  one  stationary 


and  one  travelling 


eiqi.r^1* 


(15) 


£2(q2 ,r  cos  wTt  -  q2,»  sin  u)Tt)ea2t 


(16) 


linear  disturbances  being  present  in  the  signal,  one  may  first  extract  information  for  the 
damping  rate  of  the  stationary  mode  from  the  signal  itself  and  for  the  damping  rate  and 
frequency  of  the  travelling  disturbance  from  the  first  time-derivative  of  the  DNS  signal  for 
q.  Subsequently,  one  may  solve  the  (2 NxNy)  x  (2 NxNy)  system  defined  by  writing 


(17) 


[°2  +  )q  +  (a^  +  cr2  +  w*  —  2aia2)eiqi,recrit  = 

<92q  3q  2  2. 

__2(r2__  +  (<72+air)q 


at  two  consecutive  times  fi  and  £2  for  q  and  £iqit7.,  where  Nx  and  Ny  are  the  number 
of  points  discretising  the  x—  and  y— spatial  directions,  respectively. 

The  accuracy  by  which  wT  and  a  are  determined  depends  on  that  by  which  the  first 
three  time-derivatives  of  q  are  calculated;  this,  in  turn,  depends  on  the  time-step  in  the 
calculation  and  the  number  of  fields  stored  in  order  for  backward  differentiation  formulae 
to  be  applied.  Since  the  time-step  is  controlled  by  CFL  considerations,  it  is  advisable  to 
store  a  reasonably  high  number  of  fields  in  order  for  high  accuracy  of  and  a  and,  in  turn, 
of  q  to  be  obtained.  The  calculations  to  be  presented  in  what  follows  have  been  performed 
using  five-point  backward  differencing  formulae  on  an  equidistant  grid  [14]. 

At  conditions  at  which  a  steady-state  solution  exists  most  two-dimensional  global  eigen- 
modes  of  the  converged  steady-state  are  heavily  damped  (an  =  0(1)  in  equation  (14)). 
Consequently,  if  the  time-integration  of  the  equations  of  motion  is  pursued  long  enough, 
only  a  handful  of  (f3  =  0) -global  eigenmodes  will  survive  and  persist  in  the  DNS  signal. 
Clearly,  it  is  the  least  damped  of  the  global  instabilities  that  will  determine  the  ultimate 
behaviour  of  the  solution.  In  determining  whether  one  integrates  the  equations  of  motion 
until  all  but  the  least-damped  of  the  eigenmodes  have  subsided  in  order  to  apply  (12)  or  (13) 
or  one  recovers  q  at  an  earlier  time  from  a  signal  in  which  a  number  of  damped  eigenmodes 
still  persist  one  should  take  into  account  the  following  factors. 

First,  the  efficiency  of  the  specific  DNS  algorithm  determines  whether  the  cost  of  inte¬ 
grating  the  equations  of  motion  until  convergence  is  acceptable  at  given  flow  parameters. 
The  cost  of  computing  u;na,  intermediate  values  of  q  and  monitoring  convergence  of 


all  these  quantities,  possibly  for  several  eigenmodes,  must  also  be  weighed  against  the 
straightforward  approach  of  pursuing  the  time-integration  in  the  DNS  until  convergence. 
However,  at  all  Reynolds  numbers  studied  in  the  prototype  flow  monitored  both  a  and  u;r  of 
individual  modes  have  converged  within  the  first  quarter  to  half  of  the  total  integration  time, 
making  further  time-integration  superfluous.  While  the  integration  time  until  convergence 
is  short  at  low  Reynolds  numbers,  on  account  of  large  damping  rates  of  the  least-damped 
linear  eigenmodes,  at  increasingly  large  Reynolds  numbers  the  magnitude  of  the  damping 
rates  becomes  increasingly  smaller  and  application  of  the  ideas  exposed  in  this  section 
becomes  increasingly  attractive  in  order  for  substantial  savings  in  computing  effort  to  be 
materialised. 


3.  RESULTS  FOR  THE  SQUARE  LID-DRIVEN  CAVITY 

An  example  flow  in  which  these  ideas  may  be  illustrated  is  the  classic  lid-driven  cavity 
[2] .  In  its  function  as  a  testbed  for  numerous  algorithms  this  flow  has  generated  a  substantial 
amount  of  information  which  is  relevant  to  the  preceding  discussion.  Calculations  for  q 
were  performed  using  a  two-dimensional  spectrally-accurate  algorithm  for  direct  numerical 
simulation  of  flow  in  nonperiodic  geometries.  The  code  is  based  on  a  real-space  eigenvalue- 
decomposition  of  the  spectral  collocation  differentiation  matrices  extending  ideas  discussed 
by  Ku  et  al.  [8]  and  uses  one  member  of  the  low-storage  second-order  accurate  time- 
integration  schemes  put  forward  by  Spalart  et  al  [17].  A  spectral  algorithm  was  chosen 
in  order  for  optimal  accuracy  to  be  obtained  on  a  low  number  of  collocation  points,  the 
latter  being  dictated  by  the  maximum  number  of  points  on  which  numerical  solution 
of  the  partial-derivative  eigenvalue  problem  is  feasible  using  current  computer  technology. 
Solutions  were  obtained  using  Jacobi  polynomials  for  the  spatial  discretisation  at  resolutions 
depending  on  the  Reynolds  number  and  ranging  from  322  to  1282  spectral  collocation 
points.  The  time-steps  at  the  different  Reynolds  numbers  were  kept  well  below  those 


dictated  by  the  CFL  condition  in  order  for  reasonable  accuracy  of  the  results  of  a  to  be 
ensured.  In  view  of  our  arguments  being  based  on  nonparallel  linear  instability  analysis 
and  the  well-known  sensitivity  of  linear  instability  analysis  results  on  the  accuracy  of  the 
basic  flow,  we  first  present  a  validation  of  both  the  basic  flow  and  the  partial  derivative 
eigenvalue  problem. 


3.1.  Validation  studies 

The  accuracy  of  the  converged  steady-state  solutions  is  first  assessed  by  comparison  with 
the  established  works  of  Ghia  et  al  [27]  and  Schreiber  and  Keller  [19].  Converged  basic 
states  have  been  calculated  at  several  Reynolds  numbers  of  which  we  present  calculations 
at  Re  —  400,1000,3200  and  4000,  the  first  three  obtained  on  322  and  the  last  on  482 
Legendre  collocation  points.  At  Re  =  400  and  1000  both  aforementioned  works  present 
results  while  at  the  higher  Reynolds  numbers  we  compare  our  calculations  individually  with 
either  work.  Interestingly,  aside  from  the  locations  and  maximum  values  of  streamfunction 
and  vorticity  in  the  primary  vortex  core,  Schreiber  and  Keller  [19]  analysed  and  presented 
their  results  in  the  form  of  a  converging  series  calculated  by  Richardson  extrapolation. 
Comparisons  are  presented  in  a  twofold  manner.  The  comparison  of  our  calculations  for 
the  location  and  maxima  in  the  stream-function  ^  and  the  vorticity  £  with  those  of  the 
reference  works  is  shown  in  Table  la;  note  that  [19]  define  (  to  have  an  opposite  sign  to 
that  of  [27]  and  the  present  work.  Although  the  overall  agreement  of  all  results  is  quite 
reasonable  marginal  differences  exist.  These  may  be  attributed  to  the  different  grids  used  in 
all  three  works,  making  an  interpolation  procedure  necessary  for  detailed  comparisons.  To 
this  end,  we  employed  a  piecewise  cubic  procedure  to  transfer  our  results  onto  the  (different) 
maxima  of  the  benchmark  calculations.  Our  interpolated  values  as  well  as  the  results  of 
[27]  and  [19]  are  presented  in  Table  lb,  where  the  individual  comparisons  demonstrate  a 
substantially  more  satisfactory  agreement  of  our  calculations  with  both  benchmark  works 


at  low  Re- values  and  especially  with  the  Richardson  extrapolated  results  of  [19]  at  the 
highest  Reynolds  number  monitored. 

It  is  well-known  from  comparisons  of  three-dimensional  DNS  results  and  one-dimensional 
Orr-Sommerfeld-based  linear  instability  analysis  that  details  of  the  steady  basic  state 
strongly  influence  the  accuracy  of  the  growth/damping  rates  of  linear  eigenmodes  [12]. 
The  remaining  differences  between  our  results  for  the  two-dimensional  steady-states  in 
the  lid-driven  cavity  and  those  of  the  benchmark  works  are  next  assessed  in  this  light, 
from  the  point  of  view  of  their  influence  on  the  global  linear  instability  analysis  results. 
Two  solutions  of  the  partial  derivative  eigenvalue  problem  (6a-6d)  for  the  lid-driven  cavity 
exist,  those  of  Ramanan  and  Homsy  [21]  (RH)  and  Ding  and  Kawahara  [5]  (DK).  These 
authors  have  presented  linear  instability  analyses  of  the  square  lid-driven  cavity  flow  which 
deliver  consistent  results  at  low  Re  but  predict  different  critical  Reynolds  number  values 
for  linear  amplification  of  three-dimensional  perturbations.  While  individual  comparisons 
are  certainly  possible,  at  high  Reynolds  numbers  neither  work  presents  results  for  the  two- 
dimensional  global  linear  instabilities  which  are  central  to  the  theme  of  the  present  paper. 
We  therefore  refrain  here  from  discussion  of  three-dimensional  linear  instability  and  the 
issue  of  a  linear  critical  Reynolds  number  and  monitor  a  low  value  of  the  Reynolds  number, 
Re  =  200,  at  which  both  RH  and  DK  present  results  at  fi  —  0. 

Table  2  shows  the  tabulated  values  of  RH,  the  graphically  reproduced  results  of  DK  and 
our  solutions  of  the  partial -derivative  eigenvalue  problem  (6a-6d).  The  overall  agreement 
of  the  previous  and  the  present  instability  analyses  is  quite  good  and  all  results  indicate  the 
experimentally  established  fact  of  stability  of  the  two-dimensional  flow  in  the  lid-driven 
cavity  at  this  Reynolds  number  [2].  Regarding  the  quality  of  the  basic  flow,  it  may  be 
inferred  from  the  results  of  Table  2  that  the  basic  states  of  both  RH  and  DK  and  the  present 


work  are  practically  identical  for  the  purposes  of  the  linear  instability  analysis  that  follows, 
at  least  at  the  Reynolds  numbers  discussed. 

3.2.  Numerical  residuals  and  ( 0  =  0)  linear  eigenmodes  in  the  square  lid-driven 

cavity 

Figs.  l-4c  show  the  convergence  histories  of  the  two-dimensional  DNS  at  several 
Reynolds  numbers,  with  the  qualitative  behaviour  of  residuals  discussed  earlier  observed. 
The  convergence  of  the  rate  of  decay  of  residuals  a,  calculated  using  (3-4),  is  shown  in 
Table  3  at  Re  —  100, 200  and  300.  Also  shown  is  the  damping  rate  vi  of  the  least- 
damped  eigenmode  having  0  =  0  as  obtained  by  linear  analysis,  based  on  the  partial- 
derivative  eigenvalue  problem  (6a-6d),  of  the  converged  steady-state  q  corresponding  to 
each  Reynolds  number.  The  excellent  agreement  between  the  two  quantities  leaves  lit¬ 
tle  room  for  doubt  that  numerical  residuals  may  be  identified  as  being  the  least-damped 
[0  =  0) -eigenmode  of  the  corresponding  converged  steady-state.  It  is  interesting  to  note 
here  that  such  an  agreement  could  not  be  obtained  when  we  followed  the  commonly-used 
procedure  to  terminate  the  steady-state  calculation  after  a  decay  of  residuals  by  an  arbitrar¬ 
ily  defined  seemingly  adequate  small  number  of  orders  of  magnitude,  say  5-6.  Such  poorly 
converged  in  time  basic  states  may  be  viewed  as  comprising  a  small  unsteady  component 
the  linear  instability  analysis  of  which  is  bound  to  deliver  erroneous  results.  Further,  it  is 
worth  mentioning  that  the  prediction  (7)  of  the  time  necessary  to  integrate  the  equations  of 
motion  until  convergence  in  time  is  in  line  with  the  results  of  Table  3  and  Fig.  1. 

A  clearly  defined  single  value  of  a  which  determines  the  behaviour  of  residuals  in  the 
entire  course  of  the  time-integration  is  a  result  of  a  two-dimensional  eigenspectrum  of  q 
in  which  the  least  damped  two-dimensional  (0  =  0)  -eigenmodes  are  stationary  and  well 
separated  in  parameter  space  from  their  more  stable  counterparts.  The  situation  becomes 
more  intricate,  but  still  amenable  to  analysis,  as  the  Reynolds  number  increases.  Qualitative 
differences  may  be  found  between  the  results  of  Figs.  1,  2  and  3,  although  all  simulations 


were  started  from  the  same  initial  condition  =  £  =  0;  we  discuss  the  differences 
between  the  results  of  Figs.  1  and  3  first.  While  in  both  sets  of  results  a  short  initial 
transient  is  followed  by  exponential  decay  of  residuals,  in  the  first  set  this  decay  pursues  at 
the  same  rate  for  almost  two  decades  while  in  the  second  two  different  rates  of  exponential 
decay  of  residuals  are  demonstrated.  Inspection  of  the  full  spectra  delivered  by  numerical 
solution  of  (6a-6d)  at  each  Reynolds  number  reveals  that  as  the  Reynolds  number  increases 
an  increasingly  larger  number  of  eigenmodes,  both  stationary  and  travelling  appear  in 
the  eigenspectrum  of  q,  having  damping  rates  approximately  equal  with  that  of  the  least 
damped  eigenmode.  As  a  consequence,  the  numerical  solution  may  initially  be  attracted 
to  a  different  than  the  least  damped  (0  =  0) -eigenmode  but  its  long-time  behaviour  will 
be  determined  by  the  latter  disturbance.  In  both  the  Re  —  500  and  the  Re  =  1000  results 
of  Fig.  3  the  damping  rates  u>i  of  the  least  and  the  next  more  stable  mode  are  presented 
as  symbols  superimposed  upon  the  curves  used  to  determine  cr.  The  results  at  Re  —  400, 
on  the  other  hand,  are  a  qualitatively  different  manifestation  of  the  same  phenomenon  of  a 
dense  eigenvalue  spectrum  in  the  neighbourhood  of  the  least  stable  (0  =  0) -eigenmode 
of  q.  Instead  of  the  solution  being  attracted  by  two  distinct  linear  eigenmodes  for  long 
integration  times,  here  the  progression  between  a  q  initially  attracted  by  the  third  least  stable 
member  of  the  eigenspectrum  of  q  to  the  final  state  (which  is  again  determined  by  the  least 
stable  eigenmode)  is  gradual,  taking  place  throughout  the  entire  convergence  history.  The 
result  is  the  barely  perceptible  deviation  from  an  exponential  decay  of  residuals  which  may 
be  seen  in  the  result  for  the  time  dependence  of  log{^}  at  this  Reynolds  number,  presented 
in  Fig.  2. 

Yet  another  qualitatively  different  behaviour  is  observed  in  the  time-signal  of  q  as  a 
consequence  of  a  further  increase  of  the  Reynolds  number.  Alongside  the  least  damped 
stationary  mode  travelling  disturbances  appear,  as  seen  in  the  results  of  Figs.  4a-4c.  In  all 


three  figures  a(t)  assumes  the  form  of  exponentially  decaying  disturbances.  However,  while 
at  the  lowest  Reynolds  number  a  clearly  identifiable  sinusoidal  perturbation  may  be  seen, 
having  u/r  «  0.97±0.01,  a  barely  perceptible  deviation  from  a  single  oscillatory  disturbance 
(wT  «  0.954  ±  0.012)  may  be  seen  at  Re  —  5000;  at  Re  =  7500  the  solution  demonstrates 
a  behaviour  which  might  be  interpreted  either  as  nonlinearity  or  as  superposition  of  two 
exponentially  decaying  linear  sinusoidal  disturbances  having  frequencies  of  u)T  =  0.933 
and  0.945.  In  order  to  analyse  these  observations  we  pursue  two  independent  paths.  First, 
we  perform  a  nonparallel  linear  instability  analysis  of  the  converged  steady  state  at  each 
Reynolds  number  and  monitor  the  least-stable  member  of  the  eigenspectrum,  which  turns 
out  to  be  a  stationary  linear  eigenmode.  Second,  we  perform  a  discrete  Fourier  transform 
(DFT)  of  the  DNS  signal  for  q  at  Re  =  2500, 5000  and  7500  and  compare  the  results  with 
the  eigenvalues  of  the  travelling  disturbances  delivered  by  the  linear  instability  analysis. 

Table  4  shows  that  a  progressive  deviation  of  the  rate  of  decay  of  the  residuals  from  the 
damping  rate  of  the  least-damped  eigenmode  occurs  as  Re  increases.  This  result  suggests 
that  as  the  Reynolds  number  increases  nonlinear  interaction  of  the  least  stable  eigenmodes 
may  cause  a  departure  of  the  numerical  solution  from  a  behaviour  predicted  by  nonparallel 
linear  theory.  The  role  that  the  least  stable  members  of  the  full  eigenvalue  spectrum  play 
in  the  dynamics  of  the  flow  may  be  inferred  from  the  results  of  Figs.  5a-  5c.  In  Fig.  5a 
we  present  the  DFT  of  the  DNS  signal  for  t/>(0.5,  0.5)  scaled  by  the  maximum  value  of  the 
spectral  density.  A  single  peak  at  2ivf  «  1,  albeit  of  somewhat  wide  support,  dominates 
over  two  much  smaller  peaks  at  2irf  =  0  and  2nf  «  2.  Shown  are  also  the  results 
of  (6a-6d)  for  u>r,  arbitrarily  placed  on  the  vertical  axis  for  readability.  An  one-to-one 
correspondence  between  the  peaks  in  the  spectrum  and  the  values  of  u>r  for  stationary  and 
travelling  linear  eigenmodes  is  clearly  identifiable.  Interestingly,  the  width  of  the  support 
of  the  peaks  is  found  to  be  associated  with  the  existence  of  more  than  one  eigenvalues  in 


the  partial-derivative  eigenvalue  problem  spectrum,  at  both  lot  «  1  and  wr«2.  The  origin 
of  the  existence  of  only  harmonics  of  the  first  travelling  eigenmode  in  the  full  eigenvalue 
spectrum  deserves  further  investigation.  It  should  be  noted  here  that  a  Krylov  subspace 
iteration  method  has  been  used  for  the  solution  of  the  partial-derivative  eigenvalue  problem, 
which  results  in  only  a  window  of  the  eigenvalue  spectrum  being  captured  at  any  single 
calculation.  The  number  of  converged  eigenvalues  recovered  increases  as  the  subspace 
dimension  increases.  However,  the  neighbourhood  of  (u;r,  a;*)  =  0  has  been  well  resolved 
in  all  results  presented  here.  A  higher  Krylov  subspace  dimension  has  been  found  to  deliver 
additional  eigenvalues  at  higher  frequencies. 

While  the  agreement  between  the  frequencies  in  the  DNS  signal  and  those  of  the  non¬ 
parallel  linear  analysis  of  the  converged  steady  state  is  evident,  the  results  of  Fig.  5a  do 
not  provide  any  information  on  the  damping  rates  u>i  of  the  disturbances  whose  frequency 
lies  at  u>r  «  1  in  relation  to  those  at  different  frequencies.  For  our  argument  that  the 
residuals  in  the  calculation  may  be  identified  as  the  least  stable  of  the  two-dimensional 
(ft  =  0)-eigenmodes  to  be  valid,  the  damping  rate  of  the  linear  disturbances  at  u>T  «  1 
must  be  lower  than  that  of  modes  with  higher  frequencies;  this  is  a  point  to  which  we  will 
return  shortly.  Qualitatively  analogous  results  are  obtained  at  Re  =  5000,  seen  in  Fig.  5b. 
Besides  the  slight  shift  towards  lower  frequencies,  the  quantitative  difference  with  the  re¬ 
sults  at  Re  =  2500  is  that  the  strength  of  the  eigenmodes  at  u;r  «  2  is  substantially  larger 
than  that  of  their  counterparts  at  Re  =  2500  in  relation  to  the  strength  of  the  respective 
modes  at  u>T  «  1;  this  is  the  origin  of  the  slight  deviation  from  a  purely  sinusoidal  behaviour 
of  the  signal  at  this  Reynolds  number,  seen  in  Fig.  4b.  Further,  additional  eigenmodes  ap¬ 
pear  alongside  the  counterparts  of  those  seen  at  Re  =  2500  at  vT  —  0  and  uT  «  1  and 
new  modes  appear  at  ur  «  3  and  cvT  «  4.  Finally,  at  Re  —  7500,  the  pattern  discovered 
at  the  lower  Reynolds  number  values  qualitatively  repeats  itself,  with  the  appearance  of 


additional  modes  at  the  same  and  new  at  high-frequencies;  furthermore  a  new  mode  which 
does  not  fit  in  the  period-doubling  scenario  discussed  also  is  present  in  the  linear  instability 
analysis  results,  which  the  DFT  reveals  to  be  too  weak  to  play  an  important  role  in  the 
dynamics  of  the  flow  at  this  Reynolds  number  value. 

We  return  to  the  question  of  damping  rates  of  the  linear  instability  modes  and  present 
in  Fig.  6  the  full  spectrum  of  eigenvalues  in  the  neighbourhood  of  u)  =  0  at  Re  = 
2500, 5000,  and  7500.  Results  of  significance  in  this  figure  are  the  following.  First,  as 
the  Reynolds  number  increases  the  flow  becomes  less  stable  to  two-dimensional  linear 
(/3  =  0)-eigenmodes.  Second,  in  all  three  Reynolds  numbers  the  least  stable  modes 
are  stationary  disturbances.  Third,  perfect  symmetry  about  u;r  —  0  may  be  observed 
in  the  results,  as  should  be  expected  from  the  ability  to  reformulate  (6a-6d)  as  a  real 
eigenvalue  problem.  Consistent  with  the  DFT  results  of  the  signal  discussed  earlier, 
the  eigenmodes  at  «  1  are  less  stable  than  their  counterparts  at  higher  frequencies. 
Comparing,  for  example,  the  Re  —  5000  eigenmodes  (u>r,u;i)  =  (0.967,-0.0158)  and 
(a>r,Ci>i)  =  (1.921,  -0.0319)  one  finds  that,  if  introduced  at  the  same  initial  amplitude  in 
the  flow,  the  second  mode  would  be  reduced  by  a  given  number  of  orders  of  magnitude  in 
amplitude  in  approximately  half  as  long  an  integration  time  as  that  required  for  the  first 
mode  to  experience  the  same  reduction  of  amplitude. 

3.3.  The  critical  Reynolds  number  of  (0  =  0)— linear  disturbances 

The  preceding  discussion  leads  to  re-examination  of  the  question  of  a  critical  Reynolds 
number  for  linear  growth  of  two-dimensional  global  instabilities  in  the  square  lid-driven 
cavity.  Consistent  with  well-established  numerical  solutions  for  the  steady-state  in  this  flow 
the  nonparallel  linear  instability  analysis  results  of  §  3.2  deliver  a  least  damped  stationary 
({3  =  0)  -eigenmode  which  has  a  damping  rate  whose  magnitude  decreases  with  increasing 
Reynolds  number.  The  dependence  of  on  Re  for  this  mode  has  been  obtained  at  several 


Reynolds  numbers  and  is  presented  by  symbols  in  Fig.  7.  Analysis  of  the  results  for  the 
damping  rate  of  the  least-damped  eigenmode  as  function  of  the  Reynolds  number  delivers 
a  curve-fit  of  the  data  by  using 


u>i  =  —109.071  Re~im8.  (18) 

The  curve  defined  by  (18)  is  also  shown  in  Fig.  7  by  a  solid  line  and  has  been  found  to 
deliver  reasonably  accurate  predictions  of  at  Re  >  1000,  where  the  calculated  data  may 
be  collapsed  onto  a  single  curve.  The  upper  bound  of  the  Reynolds-number  range  in  which 
(18)  may  be  used  with  confidence  to  predict  the  rate  of  decay  of  residuals  and  the  associated 
time  of  integration  of  the  equations  of  motion  until  a  steady-state  solution  is  reached  must 
be  Re  «  104,  a  value  below  which  a  multitude  of  two-dimensional  numerical  solutions 
have  demonstrated  the  existence  of  converged  two-dimensional  states.  In  the  framework  of 
the  current  nonparallel  linear  instability  analysis  this  should  manifest  itself  by  q  losing  its 
stability  in  a  linear  framework  to  amplified  two-dimensional  perturbations  having  f3  =  0. 
However,  as  has  been  mentioned,  the  existence  of  a  converged  steady-state  solution  is 
synonymous  with  all  global  eigenmodes  of  the  flow  being  stable.  Another  possibility  is 
that  the  nonlinear  interaction  of  two-dimensional  global  neutrally-stable  disturbances  as 
the  Reynolds  number  increases  may  be  held  responsible  for  the  observed  inability  to  obtain 
a  converged  steady  state  solution.  However,  from  (18)  it  follows  that  Ui  <  0,  VJf?e  and 
the  flow  remains  stable  to  all  two-dimensional  (/ 3  =  0)— eigenmodes.  Two  aspects  of  this 
prediction  should  be  stressed  here.  First,  (18)  is  a  curve-fit,  at  best  valid  up  to  the  highest 
Reynolds  number  used  to  produce  it,  Re  =  7500.  Second,  the  filling-up  of  the  eigenspec- 
trum  and  the  associated  nonlinear  interaction  of  some  of  the  least  stable  eigenmodes  as  the 
Reynolds  number  increases,  causes  a  systematic  departure  of  the  numerical  solution  for  q 
from  one  determined  by  a  single  eigenmode  of  the  nonparallel  linear  instability  theory,  as 
already  shown  in  the  results  of  Table  4.  On  the  other  hand,  the  trend  predicted  by  (18)  is 


correct,  namely  that  the  damping  rates  of  two-dimensional  global  linear  instabilities  as  Re 
increases  are  exponentially  small  in  magnitude.  As  such,  an  increasingly  large  number  of 
global  modes  may  be  considered  neutrally  stable  at  large  Reynolds  numbers;  it  is  there¬ 
fore  likely  that  the  second  scenario,  namely  nonlinear  interaction  of  near  neutrally  stable 
two-dimensional  global  flow  eigenmodes  is  responsible  for  the  observed  loss  of  ability  to 
obtain  a  steady-state  solution  of  the  equations  of  motion  at  Re  >  104  (f.e.  [31]). 

3.4.  The  spatial  distribution  of  the  least-damped  global  eigenmode 

A  closer  look  at  the  spatial  distribution  of  the  least-damped  of  the  global  linear  eigen¬ 
modes,  which  we  have  identified  as  the  residual  of  the  DNS  calculation  at  low  and  moderate 
Reynolds  numbers,  is  now  in  order.  Before  presenting  results  of  the  eigenvalue  problem  at  a 
single  Reynolds  number  value,  Re  =  1000,  we  remark  that  the  converged  two-dimensional 
flow  q  in  the  lid-driven  cavity  has  two  velocity  components  u  and  v  both  of  which  lie  on 
the  plane  H.  Consequently,  there  is  no  preferential  direction  along  the  £ -spatial  coordinate. 
Mathematically,  this  physical  fact  is  expressed  by  the  ability  to  write  (6a-6d)  as  an  eigen¬ 
value  problem  with  real  coefficients  which  admits  either  real  or  complex-conjugate  pairs  of 
solutions.  The  former  are  identified  as  the  stationary  and  the  latter  as  the  travelling  modes. 
In  the  case  of  flow  at  Re  =  1000  the  least  damped  eigenmode  is  a  stationary  disturbance 
whose  eigenvalue  is  cj  =  (0.0,  —0.068).  Fig.  8  shows  the  scaled  disturbance  eigenvector 
as  contour-lines,  the  labelling  of  which  is  to  be  found  in  Table  5. 

The  most  striking  feature  of  the  results  for  u  and  v  is  that  these  disturbance  velocity 
components  qualitatively  correspond  to  small-perturbations  of  their  basic  flow  counterparts 
u  and  u,  respectively.  The  predominant  motion  set  up  by  the  disturbance  flowfield  follows 
that  of  the  basic  flow  with  flow  being  driven  along  the  positive  x— axis  by  u  in  the  top 
half  of  the  cavity  and  along  the  opposite  direction  in  the  lower  half.  Analogously,  v 
divides  the  cavity  in  two  parts;  one  in  which  flow  is  along  the  positive  and  one  along  the 


negative  y-direction;  this  result  may  be  visualised  in  Fig.  8.  Another  interesting  result 
clearly  visible  in  both  the  u  and  the  v  velocity  components  is  the  existence  of  local  vortical 
motion  in  the  neighbourhood  of  all  four  comers  of  the  cavity.  As  has  been  mentioned, 
the  three-dimensionality  of  physical  space  is  respected  in  the  framework  of  the  nonparallel 
linear  instability  analysis,  with  the  existence  of  a  steady-state  solution  q  being  associated 
with  a  damped  eigenmode  which  may  possess  a  nonzero  spanwise  disturbance  velocity 
component  w\  that  of  flow  at  Re  =  1000  is  shown  in  Fig.  8.  It  should  be  mentioned 
here  that  the  spatial  structure  of  w  of  the  present  two-dimensional  (f 3  =  0) -eigenmode  is 
qualitatively  reminiscent  of  that  of  its  three-dimensional  (/3  ^  0)— counterparts  discussed 
by  Ding  &  Kawahara  [5]  at  approximately  the  same  Reynolds  number  value.  Finally,  the 
disturbance  pressure  p  also  inherits  the  qualitative  characteristics  of  the  basic  flow,  with 
remarkable  analogies  existing  between  the  basic  flow  pressure  [2]  and  the  present  linear 
instability  analysis  result.  The  flowfield  pattern  set  up  by  combination  of  the  u  and  v 
disturbance  velocity  components  is  visualised  in  Fig.  9. 

Of  the  multitude  of  results  on  spurious  and  transient  solutions  of  the  two-dimensional 
incompressible  continuity  and  Navier-Stokes  equations,  we  recall  here  those  of  Schreiber 
and  Keller  [20]  and  E  and  Liu  [30],  respectively.  The  former  investigators  demonstrated  that 
a  steady-state  formulation  of  the  governing  equations  may,  on  occasion,  deliver  unphysical 
results  which  persist  when  using  different  spatial  discretisation  schemes  but  not  at  different 
resolutions.  While  speculative,  on  account  of  the  impossibility  of  direct  comparisons  of 
results  of  steady  and  unsteady  calculations,  a  plausible  argument  may  be  put  forward  here 
that  the  combination  of  the  particular  initial  conditions  and  resolution  used  in  the  steady 
calculations  result  in  the  gross  features  of  one  or  more  of  the  global  linear  eigenmodes 
being  well  resolved  at  that  resolution.  Since  the  time-dependence  which  would  ensure 
the  exponential  decay  of  the  eigenmode  at  the  Reynolds  number  values  studied  in  [20]  is 


absent  from  the  scheme  integrating  the  steady  equations  of  motion,  the  iteration  merely 
serves  to  converge  to  a  solution  which  is  close  to  the  initial  attractor.  A  higher  resolution  is 
no  guarantee  of  preservation  of  the  same  initial  condition,  resulting  in  the  irreproducibility 
of  the  behaviour  at  the  lower  resolution  and  the  characterisation  of  the  latter  as  spurious. 
E  and  Liu  [30],  on  the  other  hand,  compare  different  schemes  for  the  spatial  discretisation 
used  within  time-accurate  solutions  and  present  results  at  given  times.  Some  of  the  results 
which  these  researchers  present  at  Re  =  104  are  at  a  rather  early  time  in  terms  of  that 
necessary  for  convergence  at  this  Reynolds  number  value  (Fig.4,  p.129,  streamlines  at 
t  =  2).  In  these  results  one  may  interpret  the  states  presented  as  superposition  of  several 
global  eigenmodes  which  can  be  present  at  an  0(1)  amplitude  in  the  solution  at  this  early 
time,  upon  the  well-known  steady-flow  streamline  pattern. 

3.5.  Obtaining  the  converged  steady-state  solution  from  non-converged  transient 

data 

The  preceding  discussion  has  demonstrated  the  association  of  residuals  in  two-dimensional 
incompressible  DNS  calculations  with  the  two-dimensional  global  linear  instability  modes 
of  the  converged  steady  state.  In  this  section  we  present  examples  of  recovery  of  steady- 
state  solutions  from  transient  DNS  data  using  this  information  and  the  algorithm  of  §  2.5. 
We  stress  that  the  applicability  of  the  algorithm  is  intimately  linked  with  the  quality  of  the 
DNS  and  the  initial  conditions  used  for  the  simulation,  since  both  determine  when,  for  what 
length  of  time  and  to  which  linear  eigenmode  the  time-accurate  solution  will  be  attracted 
in  the  course  of  the  time-integration.  Here  we  present  a  discussion  of  some  parameters 
which  affect  the  results  returned  by  the  algorithm  in  a  few  Reynolds  number  cases  of  those 
on  which  the  algorithm  was  validated. 

Results  at  Re  —  100,400  and  1000  are  shown  in  Table  6;  at  each  Reynolds  number 
we  have  performed  three  sets  of  calculations,  two  direct  numerical  simulations  and  one 
solution  of  the  partial-derivative  eigenvalue  problem.  Both  DNS  start  from  the  initial 


condition  ij)  =  £  =  0  for  the  flow  streamfunction  and  vorticity,  respectively.  On  the  one 
hand,  the  converged  ’exact’  steady  state  q  has  been  calculated  by  marching  the  equations 
of  motion  until  such  a  time  t  that  the  residuals  were  reduced  to  machine-roundoff  level, 
using  64-bit  arithmetic  and  monitoring  convergence  along  the  lines  discussed  in  §3.1.  On 
the  other  hand,  we  have  run  another  DNS  but  marched  the  equations  of  motion  until  such 
a  time  t  was  reached  at  which  a  linear  regime  was  identified  by  the  convergence  in  time 
of  wT  (when  applicable)  and  a.  The  time-marching  was  then  interrupted  and  either  (12) 
or  (13)  was  solved  for  the  respective  ’estimated’  steady-state  solution  q.  Finally,  the 
partial-derivative  eigenvalue  problem  (6a-6d)  was  solved  for  two-dimensional  disturbances 
(0  =  0)  developing  upon  q  and  the  eigenvalue  spectrum  pertaining  to  the  flow  at  each 
Reynolds  number  was  recovered.  The  results  were  compared  both  in  terms  of  the  magnitude 
of  the  relative  discrepancy  of  the  two  DNS-obtained  solutions  Aq  =  |(q  -  q)/q|  and  by 
monitoring  the  difference  between  a  in  the  second  set  of  DNS  and  u\.  Table  6  shows  the 
resolutions  and  time-steps  used  in  several  simulations,  the  time  t  at  which  a  converged 
steady-state  solution  (-0,  ()  was  obtained  by  DNS  and  the  time  t  at  which  the  damping 
rate  of  residuals  converged  to  within  a  predefined  tolerance  of  relative  discrepancy  10”6 
between  successive  values  of  a  and  the  results  for  x[)  were  calculated.  The  value  of  a  as 
well  as  the  relative  discrepancy  Aifi  =  |  —  -0) |  between  the  estimated  and  the  exact 

steady-states  is  also  shown;  the  level  at  which  the  eigenmode  being  damped  is  present  in 
the  transient  solution  at  time  i  may  be  inferred  from  A 

The  most  significant  result  of  this  table  is  the  ratio  t/i.  The  case  Re  —  100  is  typical 
of  one  in  which  the  least-stable  eigenmode  determines  the  transient  behaviour  of  the  DNS 
throughout  most  of  the  time-integration  process.  With  the  results  for  a  converging  quite 
quickly,  the  desired  converged  steady-state  may  be  obtained  at  a  time  between  a  quarter  at 
the  coarsest  and  a  fifth  at  the  finest  resolution  of  the  time  required  by  the  time-marching 


algorithm  for  the  residuals  to  be  eliminated.  The  result  for  a  is  only  marginally  affected 
by  resolution  and  time- step;  the  precise  times  at  which  a  converges  are  affected  by  a  small 
amount  when  refining  the  grid,  with  the  finest  resolution  results  converging  earlier.  In  all 
cases  use  of  the  algorithm  of  §  2.5  results  in  substantial  savings  compared  with  the  otherwise 
necessary  computing  effort.  The  spatial  distribution  of  the  difference  A0,  obtained  using 
48  collocation  points  to  discretise  each  spatial  direction,  is  shown  in  Fig.  10;  aside  from  the 
level  of  A0  it  is  interesting  to  notice  that  the  discrepancy  between  the  two  solutions  attains 
its  maximum  values  in  the  centre  of  the  cavity  and  neither  the  singularity  of  the  boundary 
conditions  nor  the  comer  vortices  are  manifested  in  this  quantity.  The  same  qualitative 
behaviour  was  shown  by  all  distributions  of  A0  at  lower  resolutions.  An  estimate  of  the 
converged  solution  0  obtained  by  application  of  (12)  at  t  —  15  may  be  found  in  Fig.  11, 
drawn  as  contours  at  the  levels  presented  by  [27].  No  cosmetic  post-processing  of  the 
results  has  been  applied,  with  values  presented  at  the  collocation  points  used.  As  it  is  to 
be  expected  by  the  results  of  Table  6  the  agreement  between  0  and  the  result  of  [27]  is 
remarkable. 

The  case  of  Re  —  400  will  be  discussed  shortly.  At  Re  =  1000,  a  converges  at 
approximately  the  same  fraction  of  total  integration  time  as  in  the  Re  =  100  results. 
However,  compared  with  the  Re  —  100  case  where  the  discrepancy  between  estimated 
and  converged  steady-states  is  three  to  four  orders  of  magnitude  smaller  compared  with 
that  shown  by  A0,  here  only  one  order  of  magnitude  difference  between  the  maxima  of 
A0  and  A0  is  shown.  Though  small,  the  discrepancy  between  0  and  0  is  much  larger 
than  roundoff  level,  implying  that  elimination  of  the  least  stable  eigenmode  from  the  time- 
dependent  signal  for  0  at  Re  =  1000  does  not  suffice  to  deliver  the  desired  0.  Another 
observation  that  may  be  made  by  comparing  the  results  of  the  Re  —  100  and  Re  =  1000 
cases  is  that  at  approximately  the  same  value  of  i/t  «  0.23,  A0  is  higher  by  about  an 


order  of  magnitude  at  Re  =  1000  compared  with  that  at  Re  =  100.  In  searching  for 
an  explanation  of  this  behaviour,  three  factors  may  be  recalled.  First,  convergence  of  a 
between  successive  values  is  a  necessary  but  not  sufficient  condition  for  the  algorithm  of 
§  2.5  to  deliver  accurate  results;  the  converged  <7  should  be  compared  with  the  corresponding 
damping  rates  in  the  least-stable  part  of  the  eigenspectrum  of  the  converged  steady-state 
q.  Second,  as  the  Reynolds  number  increases  the  damping  rates  of  all  global  eigenmodes 
decrease,  suggesting  that  increasingly  longer  integration  times  are  necessary  in  the  case 
of  a  higher  Reynolds  number  in  order  for  the  residuals  to  subside  to  the  same  level  as 
in  a  lower  Reynolds  number  case.  Third,  the  separation  of  the  eigenvalues  in  the  global 
spectrum  plays  a  significant  role  in  attracting  the  transient  solution,  as  will  be  discussed  with 
reference  to  the  Re  =  400  results.  A  distinction  must  be  made  between  the  early  and  the 
late  stages  of  the  transient  behaviour  of  the  DNS  solution.  In  the  latter  it  is  the  least-damped 
eigenmode  which  must  eventually  be  damped  in  order  for  a  steady-state  to  be  obtained. 
During  the  early  stages  of  the  simulation,  on  the  other  hand,  an  arbitrary  initial  condition 
may  need  a  large  number  of  damped  global  eigenmodes  in  order  to  be  reconstructed.  It  is, 
therefore,  conceivable  that  at  the  early  stages  of  the  simulation  a  number  of  eigenmodes 
other  than  the  least-damped  one  are  present  in  the  transient  solution.  However,  as  time 
progresses,  increasingly  more  of  these  additional  eigenmodes  subside  on  account  of  their 
large  damping  rates,  to  the  effect  that  only  the  least-damped  mode  remains  to  determine 
the  behaviour  of  the  residual.  In  other  words,  as  time  progresses,  equation  (14)  reduces  to 
(8)  and  the  theory  of  §  2.5  focussing  on  a  single  damped  eigenmode  is  applicable. 

This  conjecture  may  easily  be  put  to  test  by  simply  permitting  the  time  integration  in 
the  second  DNS  to  proceed  beyond  t  while  monitoring  on  the  one  hand  a  against 
and  on  the  other  hand  in  the  process;  the  results  may  be  found  in  Table  7.  Again 
we  discuss  the  results  at  Re  =  100  and  1000  first.  At  both  the  lower  and  the  higher 


Reynolds  number  further  integration  of  the  equations  of  motion  in  time  results  in  all  but 
the  least-stable  eigenmode  being  eliminated  from  the  signal,  as  clearly  demonstrated  by 
the  progressive  agreement  between  the  damping  rate  of  residuals  a  and  the  damping  rate 
u) i  of  the  least  stable  (/3  —  0)-global  flow  eigenmode.  Consistent  with  this  result  is  the 
increasingly  improved  accuracy  by  which  the  algorithm  of  §  2.5  returns  the  estimate  of  the 
converged  steady  state,  as  shown  by  the  minimum  and  maximum  values  of  A  ip  also  cited. 
Interestingly,  $  may  be  recovered  at  the  same  low  level  of  discrepancy  in  the  two  Reynolds 
number  cases,  f.e.  0(  10~8)  at  Re  =  100,  t/t  =  .25  and  Re  =  1000,  i/t  =  .35,  although 
the  agreement  of  a  with  w i  in  the  Re  —  100  is  about  an  order  of  magnitude  better  than  that 
in  the  Re  —  1000  case. 

At  Re  =  400,  as  has  already  been  shown  in  the  results  of  Fig.  2,  during  the  entire 
course  of  the  time  integration,  the  DNS  does  not  lock  on  any  of  the  eigenmodes  of  the 
flow.  On  account  of  the  dense  eigenspectrum  in  the  neighbourhood  of  the  least  stable 
eigenmode,  the  integration  is  first  attracted  by  the  third  least  stable  eigenmode  from  which 
it  gradually  departs  as  time  progresses,  approaching  both  the  first  and  the  second  least 
stable  eigenmodes.  In  the  meantime  continuation  of  the  time-integration  results  in  the 
converged  steady  state  being  reached  before  the  linear  eigenmodes  come  into  play  in  a 
manner  analogous  to  that  observed  at  Re  =  100  and  1000.  While  the  results  at  Re  =  100 
and  1000  suggest  that  permitting  the  time-integration  to  proceed  beyond  the  highest  value 
of  t/t  =  .35  would  yield  even  better  agreement  between  the  results  for  the  estimated  and 
the  true  steady-state  solution,  the  analogous  discrepancy  at  Re  =  400  decays  much  slower. 
Rather  than  weakening  the  findings  at  the  other  two  representative  Reynolds  numbers, 
the  Re  =  400  result  is  a  strong  demonstration  that  it  is  the  precise  details  of  the  least 
stable  part  of  the  flow  global  eigenspectrum  which  govern  the  behaviour  of  the  DNS;  if  the 
least-stable  eigenmodes  are  clearly  separated  in  the  spectrum,  as  the  case  is  at  Re  =  100 


and  Re  =  1000,  it  is  likely  that  a  single  eigenmode  will  attract  the  numerical  solution  in 
the  course  of  the  time  integration,  otherwise  a  dynamic  behaviour  of  the  DNS  which  is 
governed  at  least  for  part  of  the  integration  time  by  nonparallel  linear  instability  theory 
may  not  be  realised.  In  other  words  it  is  not  largeness  of  Re  but  the  details  of  the  global 
eigenspectrum  which  determine  whether  nonparallel  linear  theory  has  a  role  to  play  in  the 
dynamic  behaviour  demonstrated  in  the  DNS. 


3.6.  Three-dimensionality  as  a  consequence  of  amplified  (0  ^  0)  two-dimensional 

linear  eigenmodes 

Finally,  we  turn  our  attention  to  the  differences  between  two-  and  three-dimensional 
simulations  on  account  of  growing  global  linear  instability  modes.  While  the  physics 
behind  the  instability  mechanisms  is  universal,  the  lid-driven  cavity  flow  example  serves 
again  as  a  demonstrator,  with  the  differences  between  two-  and  three-dimensional  numerical 
simulation  results  in  this  flow  being  well  established  [8].  Here  we  call  upon  the  global 
linear  instability  theory  to  discuss  their  straightforward  explanation. 

It  is  possible  that  while  the  (f 3  =  0)— eigenmodes  at  a  certain  Reynolds  number  are 
damped  there  exist  unstable  0^0  global  flow  eigenmodes.  Indeed,  Ding  and  Kawahara 
[5]  have  shown  that  at  Re  =  950  the  flow  is  unstable  to  modes  having  0  €  \0i,0k]  with 
01  —  2i t/Lh  «  6.6  and  0h  ~  2tt/Li  &  8.3,  while  the  domain  of  unstable  wavenumbers 
systematically  broadens  in  both  directions  on  the 0- axis  as  the  Reynolds  number  increases. 
There  exist  two  possibilities  of  introduction  of  three-dimensionality  by  means  of  DNS, 
either  by  considering  spanwise  periodicity  (pDNS)  or  by  taking  an  aperiodic  spanwise 
domain  bounded  by  solid  walls  (aDNS).  In  the  case  of  pDNS  the  integration  domain 
in  the  spanwise  direction  is  defined  through  discrete  integer  multiples  of  a  fundamental 
wavenumber  0O  such  that  Lz  =  27r/(n/3o),  n  —  1, 2,  •  •  *,  while  in  aDNS  Lz  is  a  continuous 
free  parameter.  We  discuss  the  two  possibilities  separately;  in  both  cases  we  restrict  the 


discussion  to  simulations  performed  under  initial  and  boundary  conditions  such  that  linear 
instability  mechanisms  alone  can  drive  nonlinearity. 

If  a  three-dimensional  pDNS  is  performed  at  Re  =  950  and  a  spanwise  length  of 
the  integration  domain  Lz  is  chosen  such  that  0o  >  8.3,  that  is  Lz  <  0.76,  neither 
/30  nor  any  of  the  harmonics  of  this  global  linear  eigenmode  can  be  amplified.  As  a 
consequence  one  may  predict,  without  performing  the  three-dimensional  simulation,  that 
the  latter  will  converge  in  time  to  the  same  steady-state  solution  to  which  a  two-dimensional 
( d/dz  =  0)  simulation  converges.  At  the  same  Reynolds  number  value,  a  choice  of 
spanwise  wavelength  Lz  €  [L^Lh]  =  [0.76, 0.95]  will  result  in  exponential  amplification 
and,  eventually,  turbulent  flow  on  account  of  the  unstable  fundamental  wavenumber  which 
is  implicitly  defined  by  a  spanwise  wavelength  within  this  range.  Finally,  if  Lz  >  0.95 
two  distinct  situations  may  be  obtained;  with  the  fundamental  wavenumber  being  stable 
(J3o  <  6.6),  Lz  may  be  taken  such  that  none  of  its  harmonics  fit  within  the  domain  of 
unstable  wavenumbers  at  this  Reynolds  number,  or  an  Lz  may  be  chosen  such  that  some 
harmonic  may  be  amplified.  While  in  the  first  case  the  two-dimensional  steady-state 
solution  will  be  obtained,  the  result  of  a  three-dimensional  simulation  in  the  second  case 
will  be  transition  to  a  turbulent  flow  state.  The  case  of  an  aDNS  may  be  perceived  as 
a  special  case  of  a  pDNS,  since  the  homogeneous  Dirichlet  conditions  imposed  on  the 
disturbance  quantities  are  a  subset  of  those  admissible  in  a  periodic  simulation.  Here  there 
exist  two  possibilities,  depending  on  whether  Lz  is  smaller  or  larger  than  Li.  In  the  first 
case  two-  and  three-dimensional  simulations  will  deliver  identical  converged  steady-state 
solutions  while  in  the  second,  which  includes  the  well-studied  case  of  a  cubic  cavity, 
transition  to  turbulence  should  be  expected  on  account  of  at  least  one  three-dimensional 
(f3  ^  0)  eigenmode  having  a  wavenumber  which  fits  into  (5  G  [6.6, 8.3]  at  this  Reynolds 


number  value. 


At  higher  Reynolds  number  values  the  situation  is  qualitatively  analogous  for  aDNS, 
with  the  dichotomy  in  wavenumbers  being  determined  by  the  highest  neutrally  stable 
wavenumber  value.  For  pDNS,  on  the  other  hand,  the  analogous  discussion  to  that  at 
Re  =  950  applies  at  Lz  <  Lt  and  Lz  E  [Lh  Lh].  There  exists  a  Reynolds  number  value, 
though,  at  which  fih  >  2$;  in  such  a  situation,  if  Lz  >  Lh  there  will  always  be  some 
harmonic  of  Po  which  will  correspond  to  an  unstable  mode  having  /3  =  n/3 o  £  [fihPh] 
which  will  be  liable  to  linear  amplification  in  the  three-dimensional  simulation  and  eventual 
departure  of  the  three-  from  the  two-dimensional  numerical  simulation  results.  The  lid- 
driven  cavity  with  its  large  body  of  numerical  results  is  but  one  example  of  demonstration 
of  this  behaviour. 


4.  CONCLUSIONS 

The  questions  which  gave  rise  to  the  present  work  may  be  answered  within  the  unifying 
framework  of  the  reasonably  novel  global  linear  instability  analysis  of  a  two-dimensional 
steady  solution  of  the  equations  of  motion.  Aided  by  the  results  of  a  numerically  well- 
studied  incompressible  flow  problem  we  were  able  to  attach  physical  significance  to  the 
transient  behaviour  of  two-dimensional  time-dependent  incompressible  direct  numerical 
simulation  results.  What  is  commonly  known  as  residual  in  the  simulation  is  either  the  least 
damped  two-dimensional  (/3  =  0) -linear  eigenmode  of  the  converged  steady  state  itself, 
or  can  be  related  to  a  small  number  of  the  least  damped  modes  of  the  full  eigenvalue  spec¬ 
trum.  As  the  Reynolds  number  increases,  all  global  two-dimensional  eigenmodes  become 
increasingly  less  damped,  until  a  parameter  value  is  reached  beyond  which  no  steady-state 
solution  exists.  The  physical  information  which  is  suppressed  in  two-dimensional  simula¬ 
tions  based  on  the  steady  formulation  of  the  equations  of  motion  concerns  the  dynamical 
behaviour  of  these  two-dimensional  linear  eigenmodes.  While  unsteadiness  should  not 
be  interpreted  as  amplification  of  the  global  linear  (/3  =  0)— eigenmodes,  on  the  simple 


grounds  of  the  absence  of  a  converged  steady-state  upon  which  the  latter  would  develop, 
the  process  leading  to  unsteadiness  is  directly  linked  with  the  diminishing  magnitude  of 
damping  rates  of  the  global  linear  modes  as  the  flow  Reynolds  number  increases,  and  the 
associated  prevalence  of  nonlinearity. 

When  a  steady-state  solution  exists,  the  insight  gained  from  the  association  of  the 
transient  behaviour  in  two-dimensional  DNS  with  the  results  of  the  nonparallel  linear 
instability  analysis  of  the  converged  steady-state  may  be  used  in  a  threefold  manner.  First, 
an  algorithm  may  be  constructed,  to  recover  the  steady-state  solution  from  transient  data 
taken  well  before  convergence,  thus  making  further  time- integration  of  the  equations  of 
motion  redundant.  The  algorithm,  whose  building  elements  were  presented  in  §2.5,  is  based 
on  identification  of  the  parameters  pertaining  to  the  linear  eigenmodes  which  determine 
the  transient  behaviour  of  the  solution,  namely  the  damping  rate  a  and  the  frequency 
of  the  least  stable  eigenmodes.  Results  shown  in  §3.5  on  the  example  problem  studied 
have  demonstrated  that  up  to  three-quarters  of  the  otherwise  necessary  computing  effort 
may  be  saved  by  application  of  the  theory  of  §2.5.  Second,  the  results  of  a  nonparallel 
linear  instability  analysis  of  the  converged  steady-state  can  be  used  as  a  quality  test  of  the 
obtained  solution,  if  the  latter  has  been  obtained  using  a  time-accurate  solution  approach. 
The  rate  of  decay  of  the  residual  which  ultimately  has  to  be  damped  in  order  for  a  converged 
steady-state  to  be  obtained  should  equal  the  damping  rate  of  the  least-stable  eigenmode, 
if  both  numbers  are  substantially  larger  than  zero  in  magnitude.  Disagreement  of  these 
two  quantities  indicates  that  the  obtained  steady-state  still  contains  an  unsteady  component 
which  must  be  eliminated  by  further  time-integration,  or  by  application  of  the  ideas  of  §2.5. 
Third,  the  time  necessary  for  the  reduction  of  residuals  to  machine-roundoff  level  may  also 
be  estimated  using  nonparallel  linear  instability  theory  and  is  inversely  proportional  to  the 
damping  rate  of  the  least  damped  linear  eigenmode.  Using  the  value  of  the  damping  rate 


obtained  by  extrapolation  of  data  at  lower  Reynolds  numbers  one  predicts  that  in  the  square 
lid-driven  cavity  at  Re  =  104  a  steady-state  solution,  if  one  exists,  may  be  obtained  after 
integrating  the  unsteady  equations  of  motion  for  time  in  excess  of  t  =  4000  as  calculated 
from  (7)  and  non-dimensionalised  with  the  lid-velocity  and  the  cavity  length. 

Well  before  the  flow  tends  to  lose  its  stability  to  two-dimensional  linear  eigenmodes, 
three-dimensional  (0  ^  0) -disturbances  may  be  amplified.  Depending  on  the  size  of  the 
observation  window  in  the  third  spatial  dimension,  this  amplification  of  three-dimensional 
global  disturbances  can  explain  the  differences  between  the  results  of  the  two-  and  three- 
dimensional  DNS.  Again,  caution  is  warranted  at  this  point  not  to  confuse  amplification 
of  the  global,  two-dimensional  instabilities  discussed  here  with  solutions  of  the  classic 
ordinary-differential-equation  based  eigenvalue  problem,  which  are  incorporated  in  those 
of  (6a-6d);  both  mechanisms  may  provide  amplification,  as  the  laminar  separation  bubble 
flow  example  has  clearly  demonstrated  [29].  Conversely,  nonparallel  linear  instability 
theory  provides  a  handle  to  probe  into  the  physics  of  the  flow  in  (three-dimensional) 
physical  space  using  two-dimensional  DNS  results,  before  resorting  to  computationally 
intensive  three-dimensional  spatial  DNS,  at  least  as  far  as  the  response  of  the  flow  to  small- 
amplitude  excitations  is  concerned.  Solution  of  the  partial-derivative  eigenvalue  problem 
not  only  answers  the  question  whether  new  physics  is  to  be  learnt  by  performing  the  three- 
dimensional  DNS  at  a  given  set  of  parameters  but  also  provides  information  on  the  physical 
mechanism  which  leads  flow  to  deviate  from  two-dimensionality. 

5.  EPILOGUE 

Based  on  the  findings  presented  we  may  extend  the  discussion,  in  the  form  of  proposed 
future  work,  to  both  one  and  three  nonperiodic  spatial  directions.  Both  an  one-dimensional 
and  a  three-dimensional  steady-state  solution  q  may  be  recovered  by  application  of  the 
ideas  discussed  herein  for  the  case  of  two  nonperiodic  spatial  directions.  In  the- case  of  an 


one-dimensional  profile  q  being  seeked  by  time-marching  the  equations  of  motion,  taking 
two  spatial  directions  as  periodic  and  resolving  the  third,  the  associated  linear  instability 
problem  to  be  solved  is  based  on  the  classic  system  of  the  one-dimensional  Orr-Sommerfeld 
and  Squire  linear  instability  equations  to  which  (6a-6d)  reduce  if  the  dependence  of  the  basic 
flow  on  one  of  the  two  resolved  spatial  directions,  say  x ,  is  neglected  such  that  this  spatial 
direction  may  be  taken  as  homogeneous  as  far  as  the  disturbance  field  is  concerned.  The 
linear  mode  associated  with  the  residuals  is  the  least  stable  member  of  the  spectrum  obtained 
at  a  =  0  =  0,  a  and  0  being  the  wavenumbers  along  the  periodic  spatial  directions,  x 
and  z.  It  is  well  appreciated  in  this  case  that  agreement  of  the  time-accurate  simulation 
results  and  those  of  the  one-dimensional  linear  instability  problem  is  a  minimum  simulation 
quality  criterion  [12,  3].  However,  given  current  hardware  capabilities,  it  is  likely  that  an 
one-dimensional  q  will  be  seeked  by  a  direct  algorithm,  rather  than  by  time-marching  the 
unsteady  equations  of  motion.  An  extension  of  the  algorithm  presented  for  the  recovery  of 
a  two-dimensional  q  is  also  possible  in  the  case  of  flow  developing  in  three  nonperiodic 
spatial  directions.  In  this  case  the  existence  of  a  steady-state  q  is  synonymous  with  stability 
of  all  eigenmodes  of  the  flow  but  current  hardware  technology  makes  the  solution  of  the 
corresponding  three-dimensional  partial  derivative  eigenvalue  problem  impractical.  On  the 
other  hand  the  ideas  presented  in  §  2.5  may  be  used  in  order  to  recover  a  three-dimensional 
steady  state  once  a  regime  of  linear  damping  of  residuals  has  been  identified.  The  discussion 
presented  may  also  be  straightforwardly  extended  in  compressible  flow,  in  the  absence  of 
discontinuities.  While  two-dimensional  numerical  simulations  for  compressible  flow  in  the 
presence  of  discontinuities  are  well  advanced,  the  corresponding  linear  instability  theory 
is  only  recently  slowly  emerging  [32].  Both  nonparallel  linear  instability  theory  in  three 
nonperiodic  spatial  directions  and  that  in  the  presence  of  discontinuities  are  worth  pursuing 


in  future  studies. 
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TABLE  la 

Location  and  value  of  the  maxima  of  the  primary  and  the  lower-left  (LL),  lower-right 
(LR)  and  upper-left  (UL)  secondary  vortices  in  the  steady  state  solution 
for  ^  and  f  at  Re  =  400, 1000,  3200  and  4000;  comparisons 
with  Ghia  et  al  [27]  and  Schreiber  and  Keller  [19]. 


Re  =  400 

Ghia  et  al  [27]  Schreiber  and  Kelier  [19]  present  results 
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i> 

-0.1139 

-0.1140 

-0.1139 
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2.29469 

2.281 
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(0.5547,0.6055) 
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LL 
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Re  =  1000 

Ghia  al  [27] 

Schreiber  and  Keller  [19] 

present  results 
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-0.117929 

-0.11603 

-0.118902 
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2.02600 

2.068251 

(®>y) 

(0.5313,0.5625) 
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-1.15465 
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LL 

LR 

[27] 

if) 
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TABLE  lb 

Comparison  of  the  interpolated  values  of  our  solutions  on  the  maxima  presented 
by  Ghia  et  al  [27]  and  Schreiber  &  Keller  [19].  An  asterisk  denotes 
Richardson-extrapolated  data  in  the  latter  work. 
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i— 1 

c 

[27] 

-0.113909 

2.29469 
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-0.4335 
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-0.113989 

2.29463 

1.47210 
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6.42406 
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[27] 
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2.06839 
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-0.36575 
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-0.11894* 

2.0677* 

2.1700 
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-0.9990 
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-0.118905  : 

2.068234 

2.3151 
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-1.0481 

Re  =  3200 
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-0.121777 

1.9612 
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-1.00607 
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-2.25511 

Re  =  4000 
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UL 
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[19] 

-0.12202* 

1.9498* 

1.1200 

-1.0670 
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-0.122026 

1.94960 

1.2411 

-1.1427 

2.9228 

-2.31944 

TABLE  2 

Comparison  of  the  least  stable  eigenmode  at  Re  =  200  against  the  results  of 
Ramanan  and  Homsy  [21]  (RH)  and  the  graphically  (digitally)  reproduced 
growth  rate  result  of  Ding  and  Kawahara  [5]  (DK). 
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RH 

Wr 
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DK 
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present  results 

Wr  0>i 
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—0.34 

-0.3183 

±0.0000 
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-0.4013 

8 
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-0.5473 

TABLE  3 

Numerical  results  for  the  rate  of  decay  of  residuals  <r  as  a  function  of  resolution 
at  different  low  Reynolds  numbers.  Also  shown  the  result  of  numerical 
solution  of  (6a-6d)  for  the  imaginary  part  of  the  eigenvalue  a>i, 
using  the  respective  converged  steady-state  as  basic  flow. 


Resolution 
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Re 
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O' 
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O' 
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-0.5404 

-0.3248 

-0.2865 

24  x  24 

-0.5407 

-0.3319 

-0.2843 

32  x  32 

-0.5409 

-0.3318 

-0.2842 

40  x  40 

-0.5409 

-0.3318 

-0.2842 

Wi 

-0.5410 

-0.3319 

-0.2845 

TABLE  4 

Numerical  results  for  the  damping  rate  of  the  least  stable  (/3  =  0)— eigenmode 
at  Re  —  2500,  5000  and  7500  and  its  discrepancy  in  percentage  terms 
from  the  rate  of  decay  of  residuals  or. 


Re 

Wi 

X  100 
kl 

2500 

-0.0253 

1.2 

5000 

-0.0112 

8.9 

7500 

-0.0093 

17.8 

TABLE  5 

Contour  levels  for  the  spatial  structure  of  the  normalised  eigenfunctions  pertaining 
to  the  least  stable  global  (/ 3  =  0— )eigenmode  at  Re  =  1000,  whose 
eigenvalue  is  u>  =  (0.0,  —0.068);  x(y)  =  x  x  10y. 


u 

level 

symbol 

V 

level 

symbol 

U) 

level 

symbol 

J5 

level 

symbol 

+/-•(-!) 

A/a 

+/-s>(-i) 

A/a 

.9 

a 

.90 

a 

+/-«(-!) 

B/b 

+/-«(-!) 

B/b 

.7 

b 

.80 

b 

+/  —  3( — 1) 

C/c 

+/-*(-!) 

C/c 

.5 

c 

.70 

c 

+/  -  «(-2) 

D/d 

+/-6(-2) 

D/d 

.3 

d 

.60 

d 

+/  —  2(— 2) 

E/e 

+/-2(-2) 

E/e 

+/-•! 

e/g 

.45 

e 

+/- l{-3) 

F/f 

+/  -  l(-3) 

F/f 

0 

f 

.30 

f 

-.2 

h 

.15 

g 

0  h 

-.15  i 

-.20  j 


TABLE  6 

Recovery  of  0  from  transient  data  at  Re  =  100,400  and  1000  as  function 
of  resolution  and  time-step  used  in  the  DNS.  x(y)  =  x  x  10 v. 


Re  =  100 

Resolution 

16  X  16 

24  x  24 

32  x  32 

At 

0.01 

0.01 

0.005 

t 

50.43 

50.42 

49.005 

t 

12.71 

12.79 

11.07 

—cr 

0.540246 

0.540214 

0.540876 

max(A^) 

5.3(— 8) 

8.8(— 7) 

4.6(— 6) 

min(AV») 

3.6(— 9) 

7.9(— 8) 

5.1(-7) 

max(A0) 

3.4(— 4) 

3.5(-4) 

1.0(— 3) 

min(A^i) 

1.5(— 5) 

1.2(— 3) 

1.6(— 2) 

Re  =  400 

Resolution 

16  X  16 

24  X  24 

32  x  32 

At 

0.01 

0.01 

0.01 

i 

117.82 

110.91 

110.27 

t 

26.43 

23.43 

23.43 

—a 

0.267091 

0.258176 

0.258280 

max(AV0 

5.2(— 5) 

3.6(— 5) 

1.3(— 5) 

min(AV>) 

7.0(— 6) 

1.9(— 5) 

7.7(— 6) 

max(A-0) 

2.8(— 3) 

1.2(— 2) 

9.5( — 3) 

min(AV>) 

8.0(— 4) 

2*1(— 3) 

2.4(— 3) 

Re  =  1000 

Resolution 

24  x  24 

32  x  32 

40  X  40 

At 

0.01 

0.01 

0.01 

i 

325.93 

323.61 

324.33 

t 

77.82 

77.53 

77.81 

—(X 

0.065808 

0.065657 

0.065336 

max(A^) 

2.9(— 6) 

9.5(— 5) 

3-l(  5) 

min(A^) 

4-9(  7) 

3.2(— 6) 

3.3(— 6) 

max(AV') 

3.7(— 4) 

3.6(— 4) 

2.5(— 4) 

min(A^) 

3.7(-4) 

5.1(— 4) 

3.3(-4) 

TABLE  7 

Recovery  of  q  at  several  Reynolds  numbers  from  transient  data  at  times  beyond  that 
at  which  converges.  The  discrepancy  between  <r  and  u>i  of  the  least  stable 
eigenmode  is  also  presented.  Re  =  100  and  400  runs  on  a  32a 
grid;  Re  =  1000  run  on  a  40a  grid.  x(y)  =  x  X  10v. 


I  x  100 

max(A  ip) 

Re  =  100 

X  100 
kil 

22.59 

4.6(— 6) 

5.1(-7) 

0.023 

25.02 

4.0(— 8) 

1.9(— 8) 

0.018 

30.90 

2.5(— 9) 

1.2(-9) 

0.011 

35.12 

1.8(— 10) 

1.0(-10) 

0.009 

Re  =  400 

21.25 

1.3(— 5) 

7.7(-6) 

19.10/11.71/6.85 

25.03 

8.4(— 6) 

3.5(— 6) 

18.88/  11.49/7.07 

30.05 

2.6(— 6) 

9.5(-7) 

18.26/  10.88/7.67 

35.61 

1.3(— 6) 

5.6(-7) 

18.04/  10.66/7.90 

Re  =  1000 

23.99 

EBE9 

3.97 

BSSfl 

Iffifl 

3.35 

30.23 

2.1(-6) 

1.2(-7) 

0.65 

1.5(-8) 

5.8( — 9) 

0.24 

FIG.  3.  Convergence  history  of  ^(0.5,  0.5)  against  time  at  Re  —  500  (upper  left)  and  its  slope  (upper 
right);  lower  left  and  right,  respectively,  the  corresponding  results  at  Re  =  1000.  In  both  cases  superimposed 
and  denoted  by  symbols  are  the  eigenvalues  of  the  two  least  stable  stationary  modes. 


FIG.  4a.  The  dependence  of  the  function  d(ln  ip1) /dt  on  time  t,  showing  the  exponential  decay  of  a  single 
travelling  mode  (ux  «  0.97  ±  0.01)  superimposed  upon  the  least  damped  exponentially  decaying  stationary 
disturbance  at  Re  =  2500. 


Re 


FIG.  7.  The  dependence  of  the  damping  rate  Wj  of  the  least  damped  two-dimensional  eigenmode  of  the 
converged  steady-state  at  a  Reynolds  number  on  Re  as  predicted  by  the  model  (18)  denoted  by  the  solid  line,  and 
as  calculated  by  numerical  solution  of  the  eigenvalue  problem  (6a-6d)  denoted  by  the  symbols. 


FIG.  8 — Continued 

Disturbance  velocity  component  w. 


o  0 


FIG.  10.  The  spatial  distribution  of  the  difference  A'0(x,  $/)  =  ^  ^  at  Re  =  100  using  Nx  —  Ny  =  48 

Jacobi  collocation  points. 


FIG.  11.  An  estimate  of  the  converged  solution  $  at  Re  =  100  obtained  by  evaluating  (12)  at  t  =  15  and 
using  Nx  —  Ny  =  48  Jacobi  collocation  points.  Iso-contours  are  drawn  at  the  levels  shown  by  Ghia  et  al  [27] 


